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! ABSTRACT 

We use a set of semi-numerical simulations based on Zel'dovich approximation, friends- 

\ of-friends algorithm and excursion set formalism to generate reionization maps of high dy- 

namic range with a range of assumptions regarding the distribution and luminosity of ionizing 
sources and the spatial distribution of sinks for the ionizing radiation. We find that ignoring 

\ the inhomogeneous spatial distribution of regions of high gas density where recombinations 

are important - as is often done in studies of this kind - can lead to misleading conclusions 

\ regarding the topology of reionization, especially if reionization occurs in the photon- starved 

regime suggested by Lya forest data. The inhomogeneous spatial distribution of recombi- 
nations significantly reduces the mean free path of ionizing photons and the typical size of 

\ coherently ionized regions. Reionization proceeds then much more as an outside-in process. 

Low-density regions far from ionizing sources become ionized before regions of high gas den- 
sity not hosting sources of ionizing radiation. The spatial distribution of sinks of ionization 

\ radiation also significantly affects shape and amplitude the power spectrum of fluctuations of 

21cm emission. The slope of the 21cm power spectrum as measured by upcoming 21cm ex- 
periments should be able to distinguish to what extent the topology of reionization proceeds 
outside-in or inside-out while the evolution of the amplitude of the power spectrum with in- 
creasing ionized mass fraction should be sensitive to the spatial distribution and the luminosity 

\ of ionizing sources. 

Keywords: intergalactic medium cosmology: theory large-scale structure of Universe. 



1 INTRODUCTION 

The reionization of neutral hydrogen is an important milestone in 
the evolution of the Universe. The epoch of reionization has re- 
ceived a major boost of attention recently due to a series of obser- 
vational advances which suggest that the process is complex and 
that the reionization of hydrogen extends over wide redshift range 
from 6 < ^ < 15 (for reviews see Choudhury & Ferrara 2006a; 
Furlanetto, Oh, & Briggs 2006). We are about to enter an exciting 
phase as planned 21cm observations are expected to settle the ques- 
tions when and how the Universe was reionized. It thus timely to 
develop more accurate and detailed analytical and numerical mod- 
els in order to extract the maximum information about the physi- 
cal processes relevant for reionization from the expected large and 
complex future data sets. 

Currently operating and upcoming low-frequency radio obser- 



vations (e.g., GMRT\ 21CMA^ MWA^ LOFAR^, SKA^) of red- 
shifted 21cm emission of neutral hydrogen should also probe the 
topology of the neutral (or ionized) regions at high redshifts. Un- 
fortunately, modelling the expected data sets is not straightforward 
because of the dauntingly wide range of physical scales involved 
and our lack of knowledge of many details of the relevant physical 
processes. 

Full numerical simulations including radiative transfer ef- 
fects are still computationally extremely challenging. Modelling 
the smallest mass haloes contributing to the ionizing emissivity 
at early epochs (with a total mass M ^ lO^M©) requires linear 
scales < O.lMpc while at the same time, the size of the simulated 
regions need to extend over 100 Mpc or more in order to probe the 
largest coherently ionized regions in the final stages of reionization. 
Despite such challenging requirements, considerable progress has 
been made in performing radiative transfer simulations of ioniza- 
tion maps of representative regions of the Universe (see e.g. Gnedin 
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2000; Ciardi, Ferrara, & White 2003; Paschos & Norman 2005; 
Iliev et al. 2006b; Iliev et al. 2006a; Iliev et al. 2007; McQuinn 
et al. 2007). Most radiative transfer simulations are, however, still 
rather limited in dynamic range and equally important also limited 
in their ability to explore the large parameter space of plausible as- 
sumptions regarding the spatial distribution and time evolution of 
the ionizing emissivity. 

This is one of the reasons why modelling the evolution of ion- 
ized regions analytically using excursion- set-like formalisms has 
become a widely used and useful tool {e.g. Furlanetto, Zaldarriaga, 
& Hernquist 2004b). Such methods are well adapted to obtain esti- 
mates of the size distribution of ionized regions for arbitrary mod- 
els of the luminosity function and time evolution of the ionizing 
emissivity. One has, however, to keep in mind that these models 
make quite drastic simplifying assumptions. The shapes of ionized 
bubbles are e.g. assumed to be spherical and the (relative) spatial 
distribution of sources and sinks of ionizing radiation are not prop- 
erly taken into account. 

Ideally, one would like to compare realistic models of the ion- 
ization state of the IGM with a large dynamic range for a wide 
range of assumptions with future observations. For this purpose a 
variety of semi-numeric formalisms have recently been proposed 
which are based on performing an excursion- set formalism on the 
initial Gaussian random field. The models predict the spatial dis- 
tribution of the (integrated) ionizing emissivity as well as the spa- 
tial distribution of ionized regions (Zahn et al. 2007; Mesinger & 
Furlanetto 2007; Geil & Wyithe 2008). They incorporate many of 
the relevant physical processes and allow the modeller to produce 
21cm maps for representative volumes of the Universe with a mod- 
est computational effort. 

These studies suggest that reionization proceeds strictly 
inside-out with dense regions ionized first and reionization slowly 
progressing into the large underdense region as time goes on 
(Furlanetto, Zaldarriaga, & Hernquist 2004b; Wyithe & Morales 
2007; Mesinger & Furlanetto 2007; McQuinn et al. 2007). This 
appears, however, in conflict with what is expected and observed 
for the post-overlap phase where the low-density regions are found 
to be highly ionized while high-density regions remain neutral be- 
cause of their high recombination rate (Miralda-Escude, Haehnelt, 
& Rees 2000; Wyithe & Loeb 2003; Choudhury & Ferrara 2005; 
Choudhury & Ferrara 2006b). These neutral regions determine 
the photon mean free path and manifest themselves as Lyman- 
limit systems in QSO absorption spectra. Based on the this low- 
redshift intuition derived from studying the intergalactic medium 
at z ~ 2 — 4 with Lya forest data, one is thus drawn to the conclu- 
sion that reionization must have proceeded - at least to some extent 
- outside-in rather than inside-out in the final stages. The obvious 
suspect for resolving this apparent contradiction is the role recom- 
binations play in these simulations. Most of these models assume a 
spatially uniform distribution of recombinations and hence do not 
take into account the self- shielding and shadowing of high-density 
regions. Furlanetto & Oh (2005) have attempted to model this by 
introducing the concept of recombination-limited bubbles in ana- 
lytic studies of the size distribution of ionized bubbles, which has 
been implemented in simulations by Lidz et al. (2008). As we will 
show in this paper it is important to realistically model the spatial 
distribution of the sinks of ionizing radiation due to recombinations 
when modelling the topology of reionization. Many of the models 
also assume that reionization occurs rather fast diminishing the rel- 
ative importance of recombinations. This appears, however, to be 
in conflict with the ionizing emissivity inferred from the opacity 
of the Lya forest in QSO absorption spectra which suggests that 



reionization occurs slowly in a photon- starved regime (Bolton & 
Haehnelt 2007). 

We will study here the effects of the inhomogeneous spatial 
distribution of recombinations on ionization maps and present a 
consistent picture of reionization combining the concepts of grow- 
ing bubbles in the pre-overlap phase with the expected presence of 
neutral clumps in the post-overlap phase. Our modelling is similar 
in spirit to other semi-numerical models of this kind and in many 
aspects we (need to) make similar approximations and simplifica- 
tions. 

In order to determine whether a high-density clump can re- 
main neutral or self-shielded against ionizing radiation, it is neces- 
sary to determine its position with respect to the nearest sources of 
ionizing photons. An important requirement for a realistic model of 
the spatial distribution of density-dependent recombinations is thus 
(i) a realistic representation of the baryon distribution and more 
importantly, (ii) the location of the sources of ionizing sources with 
respect to the density field. Note that we will here concentrate on a 
qualitative understanding of the physical effects of a spatially inho- 
mogeneous distribution of recombinations on the topology of reion- 
ization. 

The paper is organized as follows: We describe our method for 
generating the ionization maps in Section 2. Section 3 discusses our 
main results for the modelling of a single source and representative 
volumes of the Universe. In Section 4 we check the consistency 
of our modelling with Lya forest data. In Section 5 we present 
predictions for the evolution of the power spectrum and probability 
distribution of 21cm emission and discuss prospects for the first 
generation low-frequency instruments LOFAR and MWA. Section 
6 contains our conclusions. Throughout the paper, we assume a 
flat Universe with cosmological parameters Q^rn = 0.26, Qa = 
0.74, Qfe/i^ = 0.022, and h = 0.73. The parameters defining the 
linear dark matter power spectrum we use are as = 0.9, rig = 1, 
dus/d In A; = (Viel, Haehnelt, & Lewis 2006). 



2 METHOD 

Our method of constructing ionization fields at a given redshift 
consists of four steps :(i) generating the dark matter density field, 
(ii) identifying the location and size of collapsed objects (haloes) 
within the simulation box, (iii) assigning photon luminosities to 
the haloes and (iv) generating maps of ionized regions from the 
spatial distribution of the ionizing emissivity. We discuss each of 
these steps in the following subsections. 

2.1 Simulating the dark matter density field 

We obtain our representations of the dark matter density distribu- 
tion using the Zel'dovich approximation. We first generate an ini- 
tial linear density field (as is routinely done in N-body simulations) 
and then displace the particles from their initial (Lagrangian) coor- 
dinates q using the relation 

x(q,2;)=q + D+(2;)Vq0(q), (1) 

where (/)(q) is the initial velocity potential and D+ (z) is the growth 
factor of linear dark matter density perturbations. 

The advantage of the Zel'dovich approximation is its much 
larger speed compared to a typical N-body simulation of compa- 
rable size. This allows us to produce ionization maps with a very 
large dynamic range at a modest computational cost. As we will 
show later (and has been shown before) the density field obtained 
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in this way is a reasonable approximation to that obtained using full 
N-body simulations, particularly at high redshifts. 

In Figure 1 we show the projected two-dimensional density 
field of a \h~^ Mpc thick slice through our fiducial simulation. 
One can clearly identify the expected filamentary structures and 
voids, though the range of overdensities achieved at small scales are 
typically less than those obtained using full simulations. A possible 
objection against the use of the Zel'dovich approximation is that it 
becomes invalid once shell-crossing occurs. Note, however, that at 
the redshifts and at scales of our interest (> IMpc) this occurs 
rarely. A more detailed comparison of the dark matter distribution 
obtained with the Zel'dovich approximation with that of N-body 
simulations is performed in Appendix A. 

Our fiducial simulation volume is a periodic box of length 
100/i~^ Mpc (comoving) containing 1000^ dark matter parti- 
cles which corresponds to a particle mass of Mpart = 7.22 x 
lO^/i~^M0. In order to check for numerical convergence, we have 
run further simulations with differing box sizes and particle num- 
bers; these are described and discussed in Appendix B. 



2.2 Identifying haloes 

The identification of haloes within the simulation box is performed 
using a standard Friends-of-friends (FoF) algorithm (Davis et al. 
1985). Usually the mass function of haloes identified using the 
FoF algorithm with a fixed linking length h ^ 0.2 (in units 
of mean inter-particle separation) is found to give an excellent 
match to the theoretical halo mass function for masses as small as 
~ 15 — 20Mpart (for a recent example, see Springel et al. 2005). 
Unfortunately, the use of the standard linking length fails when ap- 
plied to the density field generated using the Zel'dovich approxi- 
mation due to the more diffuse matter distribution in high density 
regions. However, if we use the FoF algorithm with an adaptive 
linking length which, depending on the redshift of interest, lies in 
the range 0.30 < h < 0.37 we get very reasonable results. Note, 
that the fact that the haloes do not have the correct density pro- 
file is not a major concern here. For our purposes it is sufficient to 
obtain the correct location and mass of the haloes with respect to 
the density field. Our method is a somewhat simpler version of al- 
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Figure 2. The halo mass function at z = 6. The points with errorbars 
show the results from our simulation; the vertical errors correspond to the 
statistical uncertainties while the horizontal errors denote the bin size. The 
solid curve is the theoretical mass function of Sheth & Tormen (2002) with 
the fitting function adopted from Jenkins et al. (2001). 



gorithms based on Lagrangian perturbation theory (Monaco et al. 
2002; Scoccimarro & Sheth 2002). The main difference is that we 
identify haloes in Eulerian space using a FoF group-finder algo- 
rithm. 

The location of haloes within a slice of the simulation box 
is shown in the right panel of Figure 1. As expected the haloes 
trace the high-density peaks of the field. The comoving number 
density of haloes per unit logarithmic mass dn/dlnM at ^ = 6 
is shown in Figure 2 by the points with error-bars. The theoretical 
mass function as predicted by Jenkins et al. (2001) is shown as the 
solid curve. The agreement is excellent over a wide mass range 
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10^ < M/{h~^MQ) < 10^^. The lower mass limit corresponds 

to - lOMpart. 

2.3 Assigning ionizing luminosities 

Observationally little is known how the ionizing luminosity varies 
with galaxy properties (Inoue, Iwata, & Deharveng 2006; Chen, 
Prochaska, & Gnedin 2007; Gnedin, Kravtsov, & Chen 2008). 
Models for reionization thus often assume that the ionizing lumi- 
nosity from galaxies scales as the halo mass with an efficiency fac- 
tor chosen such that the integrated ionizing emissivity is sufficient 
to complete reionization. 

We do the same and and assume that the number of ionizing 
photons contributed by a halo of mass M is given by 



(n^(x))i? > (ni/(x))i?(l + Nr, 



(7) 



M 

N^(M)=Nion , 



(2) 



where ttih is the hydrogen mass and A^ion is a dimensionless con- 
stant. The significance of A^ion can be understood by estimating the 
globally averaged comoving photon density 



dM^iV,(M), 



(3) 



which can be written in terms of the fraction of mass in collapsed 
objects 



/coll = f dM-^M, 



riH 



(4) 



(5) 



where uh is the comoving hydrogen density. Nion is the num- 
ber of photons entering the IGM per baryon in collapsed objects 
(Wyithe & Loeb 2007). It is determined by a combination of star- 
forming efficiency within the halo, number of photons produced 
per unit stellar mass and the photon escape fraction. Note that 
the helium weight fraction Yro could equally well be absorbed 
into the definition of Mon; in that case it would be equivalent 
to the parameter ( used by Furlanetto, Zaldarriaga, & Hernquist 
(2004b) and Mesinger & Furlanetto (2007). The analysis presented 
in this paper is applicable for any functional form of N^{M). 
For example, one can include QSOs in the analysis by simply as- 
suming that they form in haloes above a given mass Mqso, i.e., 
A^7,QSo(M) = Mon,Qso(M)e(MQso/M), where 



e(x) = 1, ifx<i, 

= 0, otherwise. 



(6) 



and Mon,Qso(M) is the number of ionizing photons produced 
within a QSO-hosting halo of mass M > Mqso- 



2.4 Generating the ionization field 

Once the location and mass of haloes are known and the functional 
form of Nj{M) is assigned, the ionization field can be generated 
using an excursion- set formalism as introduced by Furlanetto, Zal- 
darriaga, & Hernquist (2004b). First we determine whether a given 
(spherical) region is able to "self-ionize". We estimate the mean 
number density of photons (n^(x))i? within a spherical region of 
radius R around a point x and compare it with the correspond- 
ing spherically-averaged hydrogen number density (ni/(x))H. The 
condition for a point x to be ionized is that 



for any R, where iVrec is the average number of recombinations per 
hydrogen atom in the IGM. For the simple model where (M) oc 
M, the above condition translates to 



(/coll (x)) J? > 



1-lHe 



(1 + iVrec), 



(8) 



which is identical to what is used in Mesinger & Furlanetto (2007). 
Points which do not satisfy the above condition are assigned a ion- 
ized fraction (5i(x) = (n-^(x))H^.^/(nif(x))i?^.^, where i^min 
is the spatial resolution of the simulation. This is important to ac- 
count for the HII regions not resolved by the resolution of the sim- 
ulations (Geil & Wyithe 2008). Note also that the effect of spatially 
uniform recombinations (i.e., the 1 + iVrec term) can be absorbed 
within the definition of A^ion- 

Before identifying ionized regions, we smooth the density 
field to a grid-size of l/i~^Mpc, corresponding to 100^ grid points 
in the box. We do this in order to smooth out the smaller scales 
which are generally comparable to the largest halo sizes where the 
Zel'dovich approximation ceases to provide a good approximation 
for the evolution of the matter distribution. 

The quantity (n^(x))i? is estimated as 



AtvR' 



(n^(x))H = f 3 



^7V^(M0e 



|x-Xi| 
R 



(9) 



where the sum is over all luminous haloes and B is defined in equa- 
tion (6). (n^(x))i? essentially measures the contribution of ioniz- 
ing photons at x arising from all the sources within a radius R 
around the point. When dealing with a small number of sources, 
the summation in the above equation can be done directly for every 
point in the simulation box. When the number of sources becomes 
large, direct summation is computationally expensive. We therefore 
convert the point source distribution into a field. The filtering is then 
done in Fourier space. The spherically-averaged hydrogen number 
density {uh (x))h is computed by assuming that the hydrogen dis- 
tribution follows the dark matter distribution and then filtering the 
density field over a scale R (Mesinger & Furlanetto 2007). 

We should mention here that our method of obtaining the ion- 
ization field follows that of Mesinger & Furlanetto (2007) with one 
notable difference. For a given R and x, we assume only the pixel at 
the centre of the sphere with radius R to be ionized when the thresh- 
old (7) is crossed while Mesinger & Furlanetto (2007) assume the 
entire filter sphere to be ionized. In this respect, our modelling is 
similar to that of Zahn et al. (2007). We have checked our method 
for isolated sources and found a good match with theoretical ex- 
pectations (to be discussed in 3.1). In the case where all (or most 
of) the identified haloes contribute to reionization, we find that the 
mass-averaged neutral fraction obtained through our method agrees 
with the theoretical value 1 — A^ion/coii/(l — ^He) to within 15 per 
cent. This difference arises because the semi-numeric schemes do 
not conserve the number of photons within overlapping ionized re- 
gions (Zahn et al. 2007). 



2.5 Implementing a more realistic inhomogeneous spatial 
distribution of sinks of ionizing radiation due to 
recombinations 

So far we have accounted for recombinations simply by multiply- 
ing the number of ionizing photons produced by a universal fac- 
tor 1 + TVrec which does not depend on location. This corresponds 
to assuming a homogeneous spatial distribution of recombinations. 
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In reality the spatial distribution of sinks of ionizing radiation due 
to recombinations will be highly inhomogeneous. Even if a given 
spherical region contains enough photons to self-ionize, the high- 
density clumps within the region will remain neutral for a longer 
period because of their high recombination rate and thus alter the 
nature of the ionization field. A simple prescription to describe 
the presence of such neutral clumps by assuming that regions with 
overdensities above a critical value (A > A^) remain neutral was 
suggested by Miralda-Escude, Haehnelt, & Rees (2000). Unfor- 
tunately for our purpose this is also not appropriate as many of 
the high-density regions are expected to harbour ionizing sources. 
Whether a region remains neutral will depend on two competing 
factors, the local density (which determines the recombination rate) 
and the proximity to ionizing sources (which determines the num- 
ber of photons available). It is thus important to include a realistic 
spatial distribution of recombinations into the formalisms for mak- 
ing ionization maps. 

As a first approximation, one can incorporate recombinations 
within the formalism by introducing a threshold condition similar 
to equation (7), i.e., 

(n^(x))H > CRaR{nH{yi))\ a"^ (10) 

where (x) is the comoving photon emissivity, an is the recom- 
bination rate at a temperature of lO^K. Note that both the number 
densities and uh are expressed in comoving units. The above 
condition, which is similar to that used by Furlanetto & Oh (2005) 
and Lidz et al. (2008), expresses the fact that for a spherical region 
of radius R to be ionized, one needs the ionizing photon emissiv- 
ity to be larger than the spherically-averaged recombination rate 
within the region. The quantity Cr is the clumping factor which 
also takes into account the fact that not all the points within the 
spherical region would contribute to the recombination rate. For 
example, Furlanetto & Oh (2005) consider that high-density points 
with A > Ai remain neutral and hence should not be counted 
while computing the recombination rate within the region. In that 
case Cr oc J^^" dAA^Py(A) is a measure of clumping factor 
provided by low-density region only, where Pv (A) is the volume- 
weighted density distribution of the IGM. 

Another possible way of modelling the recombinations in 
high-density regions is to use a self-shielding criterion. In order to 
be ionized, a given point should satisfy the condition that it cannot 
remain self- shielded, i.e., 

[n^i{yi)a-\L[yi)a]cjH < 1, (11) 

where nm (x) is the comoving number density of neutral hydrogen 
at the given point, L(x) is the comoving size of the of the absorber 
and (Th is the hydrogen photoionization cross section. In order to 
estimate the HI density for highly ionized regions, we use the pho- 
toionization equilibrium condition: nm — {aR/r)nHCi~^, where 
the photoionization rate is F = n^a~^Amfpcrif . 

Estimating F thus requires the knowledge of the emissivity 
and local mean free path Amfp. For a given filtering scale R, we 
equate it to the mean free path, i.e., Amfp = R (Lidz et al. 2008). 
Sources within a distance R then contribute to the emissivity. If 
we assume that the fluctuations in the emissivity are negligible for 
scales smaller than the mean free path, we can write 

= (^7(x)>H, (12) 
where Nj{Mi) is the photon production rate within the halo with 



mass Mi and the summation is over all haloes. The photoionization 
rate is then given by 

r(x) = {M-)Ua-'Ran cc ^ (^^) • (O) 

The above equation expresses the fact that photons travel an aver- 
age distance of R from the source before being absorbed and the 
ionizing flux at the point of the absorber is diluted by a factor R~^. 
Also implicit is the assumption that no photons are lost to recombi- 
nation within the region except those in the central cell which may 
lead to slight underestimate of the extent of self- shielded regions. 
On first sight it may appear from the above equation that sources 
which are within distances much shorter R are not properly taken 
into account (the flux from such sources would be less diluted than 
implied by tht R~^ factor). However, one has to keep in mind that 
the procedure is repeated for different values of R. Sources that 
are closer to the point will thus be taken into account for a smaller 
value of R. 

With the above approximations, one can write a new condition 
for a point to be ionized, which is 

(n^(x))H > — (14) 

There still remains the issue that the present formalism for identi- 
fying ionized regions is based on the cumulative number of pho- 
tons n^, while balancing the recombination requires the instanta- 
neous rate of photon production (as seen in the previous equa- 
tion). Any detailed model for the evolution of the ionizing emissiv- 
ity would predict both these quantities self-consistently. We would, 
however, like to incorporate recombinations here without entering 
into the complexities of the reionization history over a wide redshift 
range. We thus integrate the above from the start of reionization so 
that the left hand side gives the integrated number of photons. The 
right hand side is significant only when recombinations are impor- 
tant and hence we can write the above relation in an approximate 
way, 

(n,(x)>fl>nH(x)-^^, (15) 

^rec^XJ It 

where ti^i(x) = ai?ni/(x)a~^ is the local recombination 
timescale and etn is the timescale over which the recombination 
term has significant contribution with tn being the Hubble time. 
Note that the parameter e, which determines the time-scale over 
which recombinations are significant, depends on the ionization 
and thermal history at a given location. We also still need to ac- 
count for the effect of an enhanced recombinations due to clump- 
ing on scales smaller than resolved by our simulations. This is often 
done in the form of a sub-grid clumping factor, which should be of 
the order (but larger) than unity (Bolton & Haehnelt 2007) and can 
be absorbed within the unknown parameter e. Short of doing the 
full radiative transfer problem we have little handle for a rigorous 
estimate of e. We will thus take it to be independent of x and study 
the results for a couple of values, namely, 0.5 and 1.0. A value of 
e = 1.0 implies that recombinations are significant over a Hub- 
ble time. These values of e where chosen in order to simulate a 
model in the "photon- starved" regime of reionization suggested by 
the Lyo; forest data (Bolton & Haehnelt 2007, section 4). Smaller 
values of e should correspond to the more rapid reionization im- 
plemented in many published numerical simulations where recom- 
binations are less important. Note further that we have absorbed 
the sub-grid clumping factor within e which is typically larger than 
unity and thus e ~ 1.0 may not be that unreasonable. 
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The only parameter which remains to be discussed is the size 
of the absorber which determines the neutral hydrogen column den- 
sity. An obvious choice for this is the local Jeans length Lj{x.) 
(Schaye 2001), which depends on the temperature and density. We 
assume here a uniform temperature of lO^K. The Jeans length 
scales then asLj oc A~^/^. Any uncertainty in the value of the 
absorber size (e.g., those arising from the geometry of the object or 
a different value of temperature) would again be absorbed within 
the unknown parameter e. The ionized cells are identified using the 
two threshold conditions (7) and (15); we find that the barrier cor- 
responding to (10) is almost always weaker than (15) and thus does 
not make much difference to the results. 

Finally, we would like to point out that, while comparing the 
number of available photons to the recombination rate at a given 
point, one should exclude the collapsed gas residing within the 
halo. However, this affects only a handful number of cells within 
the box. The reason is not difficult to understand - for cells where 
the collapsed fraction is higher than, say 5-10 per cent (depend- 
ing on the exact value of A^ion being used), the cell usually pro- 
duces enough photons to ionize itself and also overcome the self- 
shielding criterion. In other words, cells with a collapsed mass frac- 
tion > 5 — 10 per cent would anyway be flagged as ionized when 
we use a filtering scale R of the order of the cell size. For cells 
with a collapsed mass fraction lower than this, it hardly makes any 
difference whether we include the halo gas into the recombination 
budget (changes of the order of a few per cent only). 

Note that our modelling probably somewhat overestimates the 
size of individual self- shielded regions. This should, however, at 
least partially be compensated by the fact that recombinations will 
occur outside of self- shielded regions and that our simulations lack 
the self- shielded regions expected to be hosted by DM haloes with 
masses below the resolution limit of our simulations. 



2.6 Other radiative transfer effects: shadowing 

There are various radiative transfer effects which have not been 
taken into account in our simplified treatment. The most important 
is the effect of "shadowing". High density clumps which are self- 
shielded from ionizing photons will not allow photon propagation 
to the other side of the source. Such shadowing effects can only 
be incorporated using some form of ray-tracing algorithm which is 
beyond the scope of the modelling here. We have studied the effect 
of shadowing for a single isolated source using a simple-minded 
ray-tracing algorithm; the details are presented in Section 3.1. 

2.7 Computational requirements 

The code used for this work has been parallelized for shared- 
memory machines using OpenMP. The simulations were run on 
COSMOS, a SGI Altix 4700 supercomputer. For our fiducial sim- 
ulation box with 1000^ particles, we used 32 processors with a to- 
tal RAM of 32 GB to store the particle positions and velocities. 
Generating the initial gaussian random field took about 10 min- 
utes, and obtaining the position and velocity data for a particular 
redshift using the Zel'dovich approximation took less than an hour. 
A substantial amount of time was required to identify the position 
of collapsed haloes using the FoF halo finder. For a single value of 
the linking length, the halo finder takes about 100 minutes to run 
using 32 processors. However, since we are using an adaptive link- 
ing length, the whole process takes much longer, about 14 hours. 
Thus, for a given redshift, generating the density and velocity fields 
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Figure 4. Evolution of the mass-averaged ionized fraction for the mod- 
els with different assumptions regarding the spatial distribution of sinks 
and sources of ionizing radiation as described in section 3.2: HR (solid 
curve), IR-0.5 (dashed curve), IR-1.0 (dot-dashed curve) and IR-HM (dot- 
ted curve). The photon production efficiency in each model is normalised 
such that xf(z = 6) = 0.8. 

alongwith the location of the haloes takes somewhere around 17 
hours. We should mention here that the FoF algorithm, which takes 
most of the time, is easily parallelizable and scales well if a larger 
number of processors is used. 

The ionization fields were generated with lower resolution. If 
we smooth the box to about 500^ grid points, the process takes 
about 40 minutes for a single set of parameters. Since we probe a 
wide range of parameter space, we usually work with a smaller 
number of grid points, say, 200^ or 100^; generating ionization 
maps takes then around a minute to complete. 



3 THE EFFECT OF SPATIALLY INHOMOGENEOUS 
RECOMBINATIONS ON THE TOPOLOGY OF 
REIONIZATION 

3.1 Test case: Stromgren sphere around a single source 
(QSO) 

First, we consider the case of a single ionizing source with the 
ionizing luminosity of a bright QSO in the most massive halo 
(M = 3.71 X 10^^ Mq) in the simulation volume as a test case. 
The ionizing luminosity and the age of the QSO are chosen such 
that the ionized region has a comoving radius of ^ 22.5/i~^Mpc 
within an otherwise completely neutral and homogeneous IGM. 
This corresponds to an ionized fraction of xf^ ^ 0.05 aver- 
aged over the whole simulation volume. Ionization maps for two- 
dimensional slices centered on the "QSO" are shown in Figure 3 
for the case with and without an inhomogeneous spatial distribu- 
tion of recombinations in the middle and left panels respectively. 
The right panel shows the same slice with the effects of shadowing 
taken into account. 

In the left panel of the figure where no (or only spatially ho- 
mogeneous recombinations have been included) the ionized region 
is spherical with radius as expected for the assumed ionizing lu- 
minosity. This confirms that our method of generating ionization 
fields is reasonably accurate for the case of spatially homogeneous 
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Figure 3. Ionization maps for a single source (QSO) with(middle panel) and without (left panel) a spatially inhomogeneous distribution of recombinations. 
The right panel includes the additional effect of shadowing. The thickness of the shce shown is l/?,~^Mpc. 



recombinations and justifies our assumption that only the central 
pixel rather than the whole filtered sphere is ionized. 

When the density dependence of recombinations are taken 
into account allowing high-density regions to stay neutral the ap- 
pearance of the ionized region is very different due to the then very 
inhomogeneous spatial distribution of sinks of ionizing radiation. 
The resulting ionized fraction also decreases from xf^ ^ 0.05 to 

^ 0.03, despite the fact that the ionizing luminosity of the 
QSO has the same value as before. This is simply due to the fact 
that a larger number of photons is needed to overcome the recom- 
binations predominantly occurring in high-density regions. More 
importantly, the shape of the ionized region is now far from spher- 
ical. The ionization fronts appear to progress into the low-density 
regions while they are halted when high-density clumps are en- 
countered (see the left panel of Figure 1 for the corresponding den- 
sity field). However, we find that there are some low-density pixels 
which lie in the shadow of a self-shielded clump but are still ion- 
ized. This unphysical "tunnelling" of photons is a limitation of our 
modelling which does not take into account shadowing effects. 

The effect of shadowing is demonstrated in the right panel of 
Figure 3. In this case we have used a simple ray-tracing algorithm 
where rays are going out from the source along all directions. For 
each point along the ray, we check whether the local photon density 
is sufficient to ionize hydrogen taking into account recombinations, 
i.e., we check whether a point can or cannot be self-shielded. The 
ray is terminated once it hits a self- shielding pixel, thus forming a 
shadow on the other side of the high-density point. The differences 
in the topology of the resulting field are obvious. The edges of the 
ionized bubble are more ragged when shadowing is included. Note, 
however, that the difference in the global ionized fraction is only 
^ 0.0003 or about ~ 1 per cent. For representative volumes of 
the Universe the effect of shadowing will be much less dramatic. 
Points lying in the shadow of a high-density clump with respect to 
one ionizing source will generally receive ionizing photons from 
sources in other directions. 



3.2 Modelling representative volumes of the Universe 

We now discuss ionization maps of representative volumes of the 
Universe where significant numbers of haloes (as opposed to a sin- 
gle source) host ionizing sources. In order to investigate the effect 
of a spatially inhomogeneous distribution of sinks and sources of 
ionizing radiation and the speed with which reionization occurs we 
consider four different models: 



• HR: The spatial distribution of recombinations is assumed to 
be homogeneous. The condition for a region to be ionized is given 
by equation (2), with Mon being chosen so as to give a defined 
global mass-averaged ionized fraction xf^. 

• IR-0.5: The spatial distribution of recombinations is assumed 
to be inhomogeneous as discussed in section 2.5. The condition for 
a region to be ionized is determined by equation (7) [which is same 
as in the HR model] and the self- shielding condition (15) with a 
value of e = 0.5. A^ion is adjusted to give the same values of x^ 
for all four models. 

• IR-LO: The same as the previous model but with e = 1.0. 
The effect of recombinations should be more prominent than in the 
previous model. 

• IR-HM: The same as model IR-1.0 except that only high 
mass (HM) haloes with M > 10^^ Mq are ionizing sources. 
This model investigates the possibility that ionizing photons within 
lower mass haloes may not be able to escape into the IGM effi- 
ciently (e.g. Gnedin 2008). These small haloes may still form stars, 
but in this model we assume that the ionizing photons are then ab- 
sorbed within the interstellar medium and hence the galaxy remains 
mostly neutral, possibly contributing significantly to the neutral hy- 
drogen budget. 

Before investigating the ionization maps generated using the 
above models, we first discuss the predicted evolution of the global 
mass-averaged ionized fraction xf^ . To obtain the evolution of xf^, 
we have calculated the collapsed mass fraction fcou using the the- 
oretical Sheth-Tormen mass function assuming a value of Mmin as 
set by our fiducial simulation box. This means that the effects of 
various feedback processes on star-formation have been ignored. 
We have then estimated the value of xf- from /coii using the rela- 
tion xf^ = Mon/coii/(l — ^He) for the HR model. For the models 
with an inhomogeneous spatial distribution of recombinations the 
above relation was modified to Xf = iVion/coll/(l - Yue) X (1 + 
etn / (tree))- The values of xf^ computed analytically in this way 
differ from those obtained using the full simulations by up to 15 per 
cent, however, the basic trends and other conclusions remain unaf- 
fected. The corresponding evolution of xf- for the four models is 
shown in Figure 4. The value of A^ion is chosen in each case such 
thatxf 0.8 atz = 6. 

For models HR, IR-0.5 and IR-1.0 (which have the same /coii 
at a given z), the evolution of xf^ is nearly identical. The growth 
of xf^ is slightly more rapid in model IR-1.0 and slightly slower in 
model HR than in model IR-0.5 , but the differences are small. For 
the same distribution of haloes, reionization progresses "faster" as 
the spatially inhomogeneous recombinations become more impor- 
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Figure 5. Ionization maps for a range of mass-averaged ionized fractions x f for the models with different assumptions regarding the spatial distribution of 
sinks and sources of ionizing radiation as described in section 3.2: HR (left-most panel), IR-0.5 (second panel), IR-1.0 (third panel) and IR-HM (fourth panel). 
The thickness of the slice shown is l/i~^Mpc. The right-most panel shows the volume-averaged ionized fraction Qi{/S) for the same models: HR (solid 
curve), IR-0.5 (dashed curve), IR-1.0 (dot-dashed curve) and IR-HM (dotted curve). 



tant. At high-z, the average recombination time is shorter than the 
Hubble time. As a result reionization is less efficient at early epochs 
when the spatially inhomogeneous recombinations are included. 
As expected the evolution of is drastically faster in model IR- 
HM, where only rare massive haloes host ionizing sources. The 
collapsed fraction in this model is significantly smaller than in the 
other models, particularly at high redshifts, and hence reionization 
is initially delayed. 

We now discuss the nature of the ionization maps for the dif- 
ferent models. Note that we have kept the halo distribution fixed 
at that corresponding to z = 6 and have varied the luminosities to 
obtain different xf' at the same redshift. In reality, however, the 



variation in is due to the evolution of the halo distribution with 
redshift. We have here chosen to keep the halo distribution fixed in 
order to focus on the effect of the different way we treat the spatial 
distribution of sinks of ionizing radiation in the different models. 

The ionization fields for different xf^ are shown in Figure 5 
with the left-most panel corresponding to model HR. The second, 
third and fourth panels describe the three models with a spatially 
inhomogeneous distribution of sinks of ionizing radiation due to re- 
combinations (IR-0.5, IR-1.0 and IR-HM, respectively). The right- 
most panel shows the volume- averaged ionized fraction Qi{/S) as 
a function of overdensity A. Including the effects of a spatially in- 
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Figure 6. Dependence of various quantities on the mass-averaged ionized fraction for the different models: HR (solid curves), IR-0.5 (dashed curves), 
IR-1.0 (dot-dashed curves) and IR-HM (dotted curves). The left panel shows the volume averaged ionized fraction x( . The dot-dot-dashed curve denotes 
xy = xf^ . The middle panel shows the number of photons produced per hydrogen atom nn in the IGM. The dot-dot-dashed curve denotes n^/nn = 
xf^ . The right panel shows the comoving mean free path Amfp of ionizing photons. 



homogeneous distribution of recombinations distinctively changes 
the topology of ionized regions at fixed ionized mass fraction. 

Let us first concentrate on the three columns of panels on the 
left of Figure 5 corresponding to models HR, IR-0.5 and IR-1.0, 
respectively. In all three models the ionizing radiation originates 
from the same dark matter haloes. The models differ only in their 
treatment of recombinations. To reach the value of , one re- 
quires higher values of A^ion in models IR-0.5 and IR-1.0 than 
in model HR as more photons are required to overcome recom- 
bination in high-density regions. When the ionized mass fraction 
is small (x^ — 0.25), the maps look very similar. At this stage 
most of the ionizing photons are ionizing the high-density struc- 
tures which host the photon sources. At the later stages of reion- 
ization {xf- > 0.5) , however, the topology of the ionized regions 
becomes very different in the three models. In model HR the topol- 
ogy of the ionized regions is significantly "smoother" than in the 
other models. The high-density regions in models IR-0.5 and IR- 
1.0 remain neutral for longer and hence a larger number of photons 
per hydrogen atom is required to reach the same x^ . Due to the 
larger number of ionizing photons per hydrogen atom the ionizing 
photons are able to reach low-density regions far away from sources 
of ionizing radiation before the average ionized mass fraction be- 
comes large. The ionization maps of model HR show much larger 
coherently ionized regions while many neutral (or partially neutral) 
clumps are embedded within the ionized regions in the models with 
a spatially inhomogeneous distribution of the sinks of ionizing radi- 
ation due to recombinations. In the very late stages of reionization, 
models IR-0.5 and IR-1.0 are nearly identical. 

The dependence of the ionization state on density is shown 
in the right-most panel; the solid, dashed and dot-dashed curves 
correspond to models HR, IR-0.5 and IR-1.0, respectively. As al- 
ready mentioned, early on (xf- < 0.5) the three models are sim- 
ilar, while they start to differ at later stages of reionization. Ini- 
tially the topology can be described as "inside-out". High density 
regions are ionized first. However, in model HR, the ionization of 
the high-density regions is fully completed before the ionization 
fronts proceed into the underdense voids, which are the last regions 
to be ionized. In the HR model reionization proceeds "inside-out" 
all the way through the reionization process. In models IR-0.5 and 
IR-1.0, on the other hand, the ionization fronts are trapped by high- 
density clumps and they therefore proceed into low-density voids 
leaving behind islands of neutral high-density gas. The topology 



is now much more complex and cannot be classified simply either 
as "inside-out" or "outside-in". Underdense regions (A < 1) are 
completely ionized by the time xf- ^ 0.75, which is expected as 
the effect of recombination is negligible within the low-density re- 
gions. For regions with A > 1, recombinations are important and 
more than one photon is required to keep the region ionized. The 
ionized fraction Qi{A) therefore decreases around A ^ 1. Higher 
overdensities A 2 are found close to the filamentary structures 
in the density field. These regions harbour small mass haloes, i.e,. 
the relatively faint ionizing sources which are able to overcome re- 
combinations to some extent and are responsible for an increase in 
the value of Qi{A) around A 2. We have verified this explic- 
itly by computing the collapsed mass fraction within such cells. 
In even higher density regions, the number of photons required to 
keep the region ionized becomes much larger than unity and can- 
not be provided by the fainter sources. Regions with overdensities 
A ^ 6 tend thus to remain neutral. The extremely high-density re- 
gions (A > 10) represent the overlapping of filaments and harbour 
the most massive/luminous sources. These regions are able to over- 
come the high recombination rates prevalent there and hence can 
remain ionized. 

The ionization maps of model IR-HM is very similar to those 
of model IR-0.5 and IR-1.0. The reversal to reionization progress- 
ing more "outside-in" occurs somewhat earlier, which is most ob- 
vious when investigating the rightmost column showing the ioniza- 
tion state as a function of density Qi{A) (dotted curves). In model 
IR-HM the ionizing sources reside in rare massive dark matter 
haloes. A significant number of high-density regions of moderate 
mass are devoid of any ionizing photon sources locally and are able 
to remain self- shielded from ionizing photons. This is very differ- 
ent from model IR-0.5 and IR-1.0 where almost all the high-density 
regions host ionizing sources and hence cannot remain completely 
neutral. Also note that the behaviour of Qi{A) for intermediate 
overdensities is somewhat different from that in models IR-0.5 and 
IR-1.0. There is no peak around A 2 in model IR-HM. Recall 
that the peak in the other models is due to the fainter sources present 
within filamentary structures. These low-mass sources are absent 
in model IR-HM and hence the corresponding peak in Qi{A) does 
not appear. 

At this point, let us briefly compare our results with other pub- 
lished results, particularly regarding the typical value of overden- 
sities which can remain self- shielded. For example, Furlanetto & 
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Oh (2005) have shown, using modeUing based on Miralda-Escude, 
Haehnelt, & Rees (2000), that an overdensity A at z = 6 would 
be self-shielded only if the local photoionization rate r-12 < 
(A/60)^/^, where r-12 is the photoionization rate in units of 
10~^^ s~^. We have explicitly verified whether this condition is 
satisfied in every self-shielded region by estimating r_i2 using 
equation (13). We find considerable fluctuations in the local value 
of r_i2 (which as expected decreases as reionization progresses 
and the mean free path rises) and there do remain regions where 
r_i2 is much lower than what is required to overcome the self- 
shielding. To give an explicit example, for the IR-1.0 scenario, we 
find that regions far away from sources have r_i2 as low as 0.002 
for xf- — 0.95 when the global mean is ^ 0.1. The range of values 
of r_i2 is typically larger than that found by Mesinger & Dijkstra 
(2008), which is probably due to the difference in the space den- 
sity of ionizing sources (the smallest haloes in the simulations of 
Mesinger & Dijkstra (2008) have a mass of ~ lO^M©, while our 
mass threshold is - 10^ M©). 

Having demonstrated that reionization should not progress in a 
simple "inside-out" manner when the inhomogeneous distribution 
of recombinations is taken into account, we now discuss various 
other quantities of interest for the different models. The dependence 
of these quantities on xf- is shown in Figure 6. 

The left panel shows the volume-averaged ionized fraction 
xY . In the HR model (solid curve) the ionized volume fraction 
does not exceed the ionized mass fraction {xj < xf^) for the 
whole range of xf^ confirming that ionization is biased towards 
high-density regions. The models with a spatially inhomogeneous 
distribution of recombinations have xY < in the early stages 
of reionization (i.e., low values of xf^), while the trend reverses 
later on. This is in line with what we discussed earlier, i.e., reion- 
ization proceeds "inside-out" at early stages while the situation is 
more complex later. As expected the reversal of trend occurs earlier 
in the IR-HM model than in the IR-0.5 and IR-1.0 models. The val- 
ues of xj are higher in the IR-HM model than in the IR-1.0 model 
for given xf^. High density regions are, on average, more neutral 
in model IR-HM, hence a larger volume has to become ionized to 
reach the same xf^ . Note that for large ionized mass fraction (say 
> 0.95) our models will increasingly underestimate the ionized 
volume fraction due to insufficient resolution. 

The middle panel shows the number of ionizing photons per 
hydrogen atom n^/nn reaching the IGM. The first point to be 
noted is that n^/riH closely follows the ionized fraction xf^ in 
model HR. Deviations arising from a moderate violation of pho- 
ton conservation of our algorithm for identifying ionized region are 
< 15 per cent. Obviously, the ratio n^/riH is higher than xf^ for 
the other models where sinks of ionizing radiation due to recombi- 
nations are included. Extra photons are required to reach the same 
ionized mass fraction. The other crucial difference between model 
HR and the other three models is that for large ionized mass frac- 
tions n^/riH flattens for model HR while it steepens when inho- 
mogeneous recombinations are included. In model HR low-density 
voids are the last regions to be ionized and hence the ionized vol- 
ume increases without significant further need for photons. The sit- 
uation is exactly opposite for the other cases where most of the 
photons are being absorbed within high-density regions (acting as 
"sinks") and hence no significant rise in xf^ is found even though 
the number of photons used up increases rapidly. 

Finally, we plot the dependence of the mean free path Amfp 
in the right panel. To calculate Amfp, we first randomly choose a 
ionized pixel and calculate the distance to a neutral pixel along a 
randomly chosen direction; this should denote the local mean free 



path for the chosen point. This Monte Carlo procedure is repeated 
for a large number of points. The global mean free path is then 
estimated in two different ways: (i) Amfp is estimated as the average 
of the different local mean free paths and (ii) Amfp is estimated as 
the median of the local mean free path distribution. In most cases, 
both methods give nearly identical estimates. The curves plotted in 
the figure are obtained using the median [method (ii)]. 

The dependence of Amfp on the ionized mass fraction is most 
easily understood in the HR model (solid curve) where it is deter- 
mined by the characteristic size of ionized regions. Amfp rises with 
x^ essentially featureless until it flattens when Amfp approaches 
the size of the simulation box. The trends for models IR-0.5 (dashed 
curve) and IR-1.0 (dot-dashed curve) are similar to that in model 
HR in the early stages of reionization {x^ < 0.4). However, as 
reionization progresses, the mean free path in models IR-0.5 and 
IR-1.0 is smaller than that in the HR model. High density clumps 
limit the propagation of ionizing photons in these models. The 
mean free path in the models with spatially inhomogeneous re- 
combination is thus not determined by the sizes of ionized regions 
when xf^ is large. It depends instead on the spatial covering factor 
of high-density peaks. Note that the mean free path in the IR-HM 
model is larger than that in the IR-0.5 and IR-1.0 models for given 
x^ . This is consistent with the fact in these models a larger volume 
has to be ionized to reach the same ionized mass fraction. Note the 
"break" in the evolution of Amfp for the models with inhomoge- 
neous recombinations. This break broadly defines the epoch when 
the mean free path starts to be limited by high-density clumps rather 
than the size of ionized regions. 

We should mention here that is likely that we have overesti- 
mated the sizes of the self-shielded absorbers because of the Hm- 
ited spatial resolution of our simulations This should lead to an 
underestimate of Amfp- The limited resolution will, however, at 
the same time, result in an underestimate of the space density of 
self-shielded regions as well as of recombination outside of self- 
shielded absorbers. This should in turn have lead to an overestimate 
of the mean free path. The two effects should thus partially cancel. 
We have examined the effect of resolution on Amfp in Appendix B 
and found that our results do not change when the resolution is im- 
proved by a factor of two. The absolute values of the mean free path 
shown in figure 6 should nevertheless be treated with some caution 
but our finding that the mean free path will evolve more slowly if 
recombinations are important should be robust. 



4 CONSISTENCY WITH LYa FOREST DATA AT z - 6 

Current observational constraints on the epoch of reionization are 
still rather limited. Studies of the Lyo; forest in QSO absorption 
spectra have taught us that reionization probably ended at around 
^ - 6 (Fan et al. 2002; Fan et al. 2004; Fan et al. 2006; c/Becker, 
Rauch, & Sargent 2007). As discussed by Miralda-Escude (2003) 
and Bolton & Haehnelt (2007), the emissivity inferred from the 
Lyo; forest data corresponds to at most a few photons per hydro- 
gen atom per Hubble time. Bolton & Haehnelt (2007) thus coined 
the term "photon- starved" to describe the regime in which reion- 
ization appears to occur. Bolton & Haehnelt (2007) measured the 
emissivity of ionizing photons to be roughly constant in comov- 
ing units in the redshift range 2 < z < 6. They pointed out that 
because of the rather low emissivity of ionizing photons reioniza- 
tion of hydrogen most likely started early and extends over a wide 
redshift range. This sits well with the rather large Thomson optical 
depth inferred from studies of the cosmic microwave background 
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(Spergel et al. 2007; Dunkley et al. 2008; Choudhury, Ferrara, & 
Gallerani 2008). Predictions of ionization maps should obviously 
be consistent with available data. Enforcing consistency with the 
Lya forest data shrinks the allowed parameter space considerably 
and we therefore discuss now how our modelling fairs in this re- 
spect. 

In Table 1 we summarize the mean-free path of ionizing pho- 
tons Amfp and the inferred photoionization rate V in our three mod- 
els for two values of the volume fraction of ionized regions (at 
z — Q) and two different assumptions for when reionization has 
started at z^e = 15 and z^e — oo, respectively. The value of 
the mean free path of ionising photons and the volume fraction of 
ionized regions x]^ at z = 6 are observationally still very uncer- 
tain. Bolton & Haehnelt (2007) estimate the mean free path to be 
<40 Mpc and infer a photoionization rate F 12 ^0.19. The vol- 
ume fraction of ionized regions has been estimated to be > 0.5 at 
2: > 6.3 from the evolution of Lya luminosity function (Kashikawa 
et al. 2006) and GRB spectrum (Totani et al. 2006), while the con- 
straints from QSO absorption line measurements at 2; < 6 are 
quoted to give 1 - > 10"^ (Fan et al. 2006). 

For our models with a spatially inhomogeneous distribution 
of recombinations the mean free path is reasonably consistent with 
the estimate of Bolton & Haehnelt (2007) if the volume fraction 
of ionized regions is large (95%). For models HR and IR-HM on 
the other hand, the estimated mean free path is consistent with the 
values in Table 1 if the volume fraction of ionized regions ai z = 6 
is low (50%). 

We have estimated the photoionization rate (in units of 10~^^ 
s~^) inferred from the photon emission rate and Amfp using 
(Bolton & Haehnelt 2007), 

where as is the spectral index of the ionizing background (which 
we assume to be 3 consistent with stellar sources of sub-solar 
metallicity). The results are shown in the two right-most columns 
in table 1. Note again that there could be inaccuracies of < 15 per 
cent arising from moderate violations of photon conservation of our 
algorithm. For models IR-0.5 and IR-1.0 we find reasonable agree- 
ment with the inferred photoionization rate for xf = 0.95. On 
the other hand, models HR and IR-HM generally tend to overpre- 
dict the photoionization rate when the assumed ionized fraction is 
large. For smaller values of the ionized mass fraction (xf = 0.5), 
these models are found to be consistent with the data. 



5 PREDICTIONS FOR 21CM OBSERVATIONS 

5.1 The effect of the spatial distribution of sinks and the 
luminosity of sources on the 21cm power spectrum. 

We have seen in Section 3.2 that the models with different assump- 
tions regarding the spatial distribution of sinks and sources of ioniz- 
ing radiation predict rather different topologies for the neutral hy- 
drogen distribution, particularly in the late stages of reionization 
(with the exception that models IR-0.5 and IR-1.0 are nearly iden- 
tical for xf^ > 0.9). We now discuss the prospects of investigating 
the effects of the spatial distribution of the sinks and sources of 
ionizing radiation with future low-frequency radio observations of 
the redshifted 21cm line. Since model IR-0.5 is qualitatively very 



Table 1. Mean free path and inferred photoionization rate for different ion- 
ized mass fractions and different redshifts for the start of reionization. 
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HR 


0.5 


15 


0.79 


0.031 


0.044 




0.95 


97 


0.96 


0.257 


0.362 


IR-0.5 


0.5 


8 


0.82 


0.016 


0.023 




0.95 


21 


1.14 


0.063 


0.089 


IR-1.0 


0.5 


5 


0.89 


0.011 


0.015 




0.95 


21 


1.40 


0.076 


0.107 


IR-HM 


0.5 


12 


0.91 


0.026 


0.036 




0.95 


56 


2.20 


0.330 


0.465 



similar to model IR-1.0, we shall not discuss it separately in this 
section. 

The 21cm brightness temperature at a given location x relative 
to the CMB can be approximated as 

Tfe(x) = nxHi(x)A(x), (17) 

where we have assumed that the spin temperature of hydrogen is 
much larger than the CMB temperature. This should be a reason- 
able assumption once a significant fraction (a few percent) of the 
volume/mass has been ionized (Scott & Rees 1990; Tozzi et al. 
2000; Ciardi & Madau 2003; Barkana & Loeb 2005; Sethi 2005; 
Furlanetto 2006; Furlanetto, Oh, & Briggs 2006; see Pritchard & 
Loeb 2008 for an extensive recent discussion of the expected evo- 
lution of the spin temperature). We have also ignored peculiar ve- 
locity effects which are small at the scales relevant here (see, e.g., 
Mesinger & Furlanetto 2007). The quantity ft ^ 22mK[(l + 
z)/7]^^'^ denotes the brightness temperature for neutral gas at mean 
density. By definition, {Tb{x.)) = fb{l-xf). 

The first quantity of interest is the power spectrum of 
temperature fluctuations which we define as T^Ali{k) = 
(T^ (k)) /27r^ . The power spectrum is plotted in Figure 7 for our 
models for a range of values of xf^ (e.g. Furlanetto, Zaldarriaga, 
& Hernquist 2004a). The panels from left to right show the power 
spectrum for models HR, IR-1.0 and IR-HM, respectively. In each 
panel the power spectrum is shown for mass-averaged ionization 
fraction xf = (solid), 0.1 (dashed), 0.3 (dot-dashed), 0.5 (dot- 
ted), 0.7 (dot-dot-dot-dashed) and 0.9 (triangles) respectively. For 
xf- = 0, the brightness temperature simply traces the DM fluctua- 
tions. 

In model HR (left panel), the amplitude of the power spec- 
trum decreases from its initial value until about xf^ = 0.3 (dot- 
dashed curve). The decrease of the fluctuation amplitude, partic- 
ularly at large scales, occurs as regions of high density become 
ionized. The decrease in amplitude is accompanied by a steepen- 
ing in slope, consistent with the findings of Lidz et al. (2007). It 
follows then a reversal in trend. The amplitude rises (particularly at 
large scales k < l.bh Mpc~^) and the slope becomes shallower. 
This is particularly evident if the power spectra for xf- —0.3 (dot- 
dashed curve) and xf- — 0.5 (dotted curve) are compared. This 
is the phase when the the ionizing radiation from collapsed objects 
ionizes the surrounding high-density regions. The growth of ion- 
ized regions boosts the large scale power and flattens the slope of 
A|i {k). The flattening of the slope continues (and the power spec- 
trum becomes practically flat) as the IGM becomes more ionized 
while the amplitude decreases at high values of xf^ . In model HR 
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Figure 7. The power spectrum of 21cm brightness temperature fluctuations. The panels from left to right show models HR, IR and IR-HM, respectively. In 
each panel, results are shown for mass-averaged ionization fraction = (solid), 0.1 (dashed), 0.3 (dot-dashed), 0.5 (dotted), 0.7 (dot-dot-dot-dashed) and 
0.9 (triangles), respectively. 
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Figure 8. The amplitude (left panel) and slope (right panel) of the power spectrum of 21cm brightness temperature fluctuations as a function of the mass- 
averaged ionized fraction xf- at wavenumber k = kp = 0.4/iMpc~^. The soUd, dashed and dot-dashed curves represent models HR, IR and IR-HM, 
respectively. 



there is nearly equal power at all scales in the late stages of reion- 
ization and the fluctuation amplitude decreases as the neutral hy- 
drogen content in the IGM decreases. 

In the early stages of reionization (x^ < 0.5) the evolution of 
the 21cm power spectrum in model IR-1.0 (middle panel) is very 
similar to that in model HR. The similarity, however, disappears 
for xf- > 0.5 when the slope of A|i (k) steepens rather than flat- 
tens. The high-density neutral regions embedded within the ionized 
regions are responsible for a considerable amount of small-scale 
power in the 21cm power spectrum. At the same time, the clumps 
limit the size of coherently ionized regions, thus keeping the large- 
scale power low. This pattern holds until the very end of reioniza- 
tion. The steepening of the slope in the later stages of reionization 
in the models with an inhomogeneous spatial distribution of re- 
combinations is a signature of the more complex topology which 
we had described in section 3.2. In the late stages reionization pro- 
ceeds much more "outside-in" than in model HR and this is clearly 
recognizable in the 21cm power spectra. 

A similar but more pronounced steepening of the slope of the 
21cm power spectrum occurs in model IR-HM, where the emis- 
sion of ionizing photons is restricted to massive haloes. Here reion- 
ization starts to proceed in a more "outside-in" fashion much ear- 
lier. The behaviour of the amplitude 21cm spectrum at large scales 
(k < l/iMpc~^) is less complicated than in the other two models; 
there is a slight dip in the power spectrum around x^ = 0.1 due to 



the ionization of the high-density regions harbouring the sources of 
ionizing radiation. Otherwise the power spectrum evolves very lit- 
tle until xf- ^0.5 and then the amplitude decreases with decreas- 
ing neutral fraction. In model IR-HM the amplitude of the 21cm 
power spectrum is generally somewhat higher than in the other two 
models. This is due to reionization being driven by relatively highly 
clustered sources in this model. 



5.2 Evolution of slope and amplitude of the 21cm power 
spectrum at scales probed by LOFAR and MWA 

In the last section we developed a feeling for how the spatial dis- 
tribution of sinks and sources of ionizing radiation influence the 
21cm power spectrum. We now discuss in more detail the possibil- 
ity to differentiate observationally between different models with 
first generation 21cm experiments like LOFAR and MWA. The typ- 
ical scales probed by these experiments correspond to wavenum- 
ber s 0.1 < /c/Mpc"^ < 1. Foreground subtraction will be a se- 
rious problem and it is not clear yet to how small and large scales 
it will be possible to determine the power spectrum with reason- 
able accuracy. We follow Lidz et al. (2007) and assume that the 
optimum scale for studying 21cm fluctuations with these instru- 
ments correspond io k ^ 0.3 — 0.5Mpc~^ and use a pivot scale of 
k — kp = 0.4/iMpc~^ in the following to be definite. This scale 
is well suited for a discrimination between our models. 
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Figure 9. The probability distribution of the dimensionless 21cm brightness temperature fluctuations. Left panel: Results for model HR . The curves with 
peaks from right to left are for mass averaged ionization fractions of xf- = 0, 0.35, 0.7, 0.9 respectively. Middle panel: Results for models IR-1.0. The 
curves with peaks from right to left represent cases xf- = 0, 0.35, 0.6, 0.9 respectively. Right panel: Results for model IR-HM. The curves with peaks from 
right to left are for xf^ = 0, 0.1, 0.45, 0.9, respectively. 



The amplitude of the power spectrum at kp and its slope rip = 
d In A|i (A:p)/d In /c as a function of the ionized mass fraction 
are shown in the left and right panels of Figure 8. The evolution 
of the different models reflects our discussion in the last section. 
The evolution of the amplitude A21 (kp) can be divided into three 
phases. An initial decrease in amplitude due to the early ionization 
of high-density regions (Wyithe & Morales 2007), is followed by a 
rise corresponding to a growth in patchiness and a final fall due to 
the elimination of neutral hydrogen. In model HR-IM (dot-dashed 
curve) where reionization is driven by rarer sources the 21cm power 
spectrum has significantly (a factor two or more) power than in the 
other two models. The amplitude of the power spectrum and its 
evolution at our pivot scale contains valuable information about the 
spatial distribution of ionizing sources. 

The amplitude of the power spectrum at our pivot scale ap- 
pears, however, not to be a good indicator of the spatial distribution 
of the sinks of ionizing radiation. Models HR and IR-1.0 models 
are very similar in this respect. As discussed earlier the spatial dis- 
tribution of sinks has instead a strong influence on the slope of the 
power spectrum Up (right panel). The evolution of the slope Up can 
be again divided into three phases. An initial rise due to the ioniza- 
tion of high-density regions followed by a fall corresponding to the 
growth of patchiness. In the third phase at xf- > 0.65 Up keeps 
on decreasing rapidly in the HR model while it increases instead in 
the models with a spatially inhomogeneous distribution of sinks of 
ionizing radiation due to recombinations. 

By measuring the power spectrum and its slope at large scales 
(k ~ 0.3 — 0.5Mpc~^) it should thus be possible to characterize 
both the spatial distribution of sources and sinks of ionizing radia- 
tion. The detectability of the 21cm fluctuations obviously depends 
on the instrument noise and the ability to subtract foreground emis- 
sion. Assuming a perfect removal of foreground emission McQuinn 
et al. (2006) find typical values of detector noise for LOFAR and 
MWA at /c - 0.4/iMpc"^ of < 5 mK^ at z = 6 for 1000 hrs 
of observation. With such noise levels, the power spectra for the 
HR and the IR-HM models should be detectable with reasonable 
confidence in the range 0.5 < x^ < 0.9 and xf^ < 0.9 ,respec- 
tively. The fluctuation amplitude in model IR-1.0 is lower than in 
the other two models discussed in this section. At their peak value 
around xf - 0.6 21 cm fluctuations should nevertheless be de- 
tectable by LOFAR and MWA even for this model. Note that the 
values quoted here should be only taken as indicative. Both the 



noise properties and the fluctuation amplitude depend on redshift 
For example in our models an ionized mass fraction of xf^ ~ 0.6 
is reached around z ^ 8 while our estimations were performed 
assuming z ^ 6. 



5.3 The PDF of the 21cm brightness distribution 

We now briefly discuss the probability distribution P{Tb/Tb) = 
P( A xhi) of the dimensionless brightness temperature (Furlanetto, 
Zaldarriaga, & Hernquist 2004b). In order to compute the distribu- 
tion, we smooth the brightness temperature Tb over scales of 10/i~^ 
Mpc, as is appropriate for the first generation 21cm experiments. 
The results are shown in Figure 9. The left panel shows the 21cm 
PDF for model HR. The curves with peaks from right to left are for 
x^ = 0, 0.35, 0.7, 0.9, respectively. We have chosen the values of 
xf^ such that they represent the characteristic points in the evolu- 
tion of the power spectrum at large scales. The curve for x^ = 
(solid) obviously represents the dark matter PDF. For xf^ = 0.35 
(dashed curve), the PDF has become significantly narrower. This 
is again due to the ionization of high-density regions and corre- 
sponds to a low- amplitude of the power spectrum. The evolution 
of the PDF in model HR is consistent with the analytical models 
of Wyithe & Morales (2007). The PDF widens subsequently with 
increasing x^ sls more regions are being ionized. The behaviour 
is similar in model IR-1.0 (middle panel) where the curves with 
peaks from right to left represent xf- = 0, 0.35, 0.6, 0.9, respec- 
tively. The only difference is a somewhat narrower width of the 
distribution than in model HR in the final stages of reionization 
(xf- ~ 0.9). This is consistent with what is expected from the evo- 
lution of the 21cm power spectra at large scales. The results for 
model IR-HM are shown in the right panel. The curves with peaks 
from right to left represent xf- — 0, 0.1, 0.45, 0.9, respectively. 
As expected, the PDF in this model is rather different from that 
in the other two models. The PDF in model IR-HM has a wider 
distribution compared to the other two models. The model predicts 
Tb/Tb > 0.5 even when the IGM is 50 per cent ionized by mass. 
Unfortunately, it is not clear whether the first generation 21 cm ex- 
periments will have enough sensitivity to constrain the shape of the 
PDF 
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5.4 Comparison with other work 

As discussed in the introduction, there has been a number of re- 
cent studies which aim at predicting the 21cm brightness distribu- 
tion. These studies range from radiative transfer simulation gener- 
ally performed by post-processing the density field of DM simula- 
tions (Ciardi, Ferrara, & White 2003; Iliev et al. 2006b; Mellema 
et al. 2006; lUev et al. 2007; McQuinn et al. 2006; McQuinn et al. 
2007; Zahn et al. 2007) to semi-numerical simulations (Mesinger 
& Furlanetto 2007; Alvarez & Abel 2007; Geil & Wyithe 2008) 
similar in spirit to the work presented here. Most of these stud- 
ies appear to agree that reionization occurs inside-out all the way 
from the start until nearly the completion of reionization. Zahn et al. 
(2007) have thereby shown that results for semi-numeric schemes 
based on collapsed mass fractions and variants of the excursion set 
formalism to identify regions which can self-ionize give very sim- 
ilar results to full radiative transfer simulations if similar assump- 
tions regarding the sources of ionizing radiation are made. Most 
similar to our work here is probably the work of McQuinn et al. 
(2007) who have studied a wide range of assumptions regarding 
the sources and sinks of ionizing radiation. When modelhng the ef- 
fects of sinks of ionizing radiation McQuinn et al. (2007) mainly 
study mini-haloes, dark matter haloes with potential wells shallow 
enough so that they can be photo-evaporated by ionizing photons. 
For these mini-haloes they find a noticeable but rather small effect 
(see Bolton & Haehnelt 2007 for a brief discussion of the role of 
mini-haloes during reionization in the photon- starved regime). Mc- 
Quinn et al. (2007), however, do not try to model recombinations 
in high-density regions in deeper potential wells which can hold on 
to photo-ionized gas in a way so that their models are likely to be 
consistent with the Lya forest data. They generally find that sinks 
of ionizing radiation and their spatial distribution have little effect 
on the topology of reionization and the power spectrum. This is 
obviously quite different from our findings. There is a number of 
differences to our modelling but the most likely reason appears to 
be the following. The emissivity used in the models of McQuinn 
et al. (2007) is rather high and reionization proceeds quickly. This 
strongly diminishes the importance of recombination compared to 
our modelling of reionization in the photon-starved regime sug- 
gested by the Lya forest data. 



6 CONCLUSIONS 

We have used here semi-numerical simulations to investigate the 
role of the spatial distribution of sinks and sources of ionizing ra- 
diation on the topology of hydrogen reionization. Our main results 
are the following. 

• The combination of Zel'dovich approximation, halo-finder 
and excursion set formalism is a powerful tool to calculate real- 
istic ionization maps with high dynamic range at a very moderate 
computational cost. 

• Enforcing consistency with the hya forest data helps to signif- 
icantly shrink the otherwise rather unconstrained parameter space 
of models of reionization. In the photon- starved regime of reion- 
ization suggested by the Lya forest data recombinations are much 
more important than in models with high ionizing emissivity where 
reionization occurs quickly. Taking into account a realistic spatially 
inhomogeneous distribution of sinks of ionizing radiation has a 
large effect on the topology of reionization in the photon- starved 
regime. 

• Initially reionization proceeds inside-out with the high-density 



regions hosting the sources of ionizing sources becoming ionized 
first. In the later stages of photon- starved reionization the sinks of 
ionizing region in our models remain neutral and reionization pro- 
ceeds deep into the underdense regions before slowly evaporating 
denser regions not hosting ionizing sources where recombinations 
are important. This reversal to a more outside-in progression in the 
late stages of reionization is more pronounced if the emission of 
ionizing radiation is restricted to massive highly-clustered and rare 
sources. 

• If the emission of ionizing radiation is restricted to rare 
sources reionization proceeds more quickly and the sizes of co- 
herently ionized regions are significantly larger. The latter results 
in an about factor two or more larger mean free path for ionizing 
photons. 

• Like other studies we find that the amplitude of the 21cm 
power spectrum and its evolution in the later stages of reioniza- 
tion is mainly sensitive to the space density of ionizing sources. 
The sensitivity to the space density of ionizing sources is, however, 
significantly increased if a realistic spatially inhomogeneous distri- 
bution of sinks of ionizing radiation is taken into account. The slope 
of the power spectrum is very sensitive to the spatial distribution of 
sinks of ionizing radiation. 

• Measurements of the amplitude and slope of the 21cm power 
spectrum at scales corresponding to A: ~ 0.3— 0.5/iMpc~^ with the 
upcoming low-frequency instruments LOFAR and MWA have ex- 
cellent prospects to reveal important information on the spatial dis- 
tribution of sinks and sources of ionizing radiation and the speed of 
reionization if the daunting tasks of accurate calibration and fore- 
ground removal are mastered successfully. The PDF of the 21cm 
brightness distribution contains important complimentary informa- 
tion. Measuring the PDF will, however, unfortunately most likely 
require higher sensitivity than can be achieved with first generation 
21cm experiments. 

Our modelling here has involved a number of significant sim- 
plifications. The spatial distribution of dark matter modelled in the 
Zel'dovich approximation was used as an proxy for the spatial dis- 
tribution of the IGM. The ionizing emissivity of sources and recom- 
bination in dense region was modelled only in an approximate inte- 
grated fashion and the dynamical effects of the ionization radiation 
on the gas were neglected. Despite the large particle number used 
in the simulations resulting in a substantial dynamic range there 
were still clear deficiencies in modelling high-density regions and 
low-mass collapsed objects/mini-haloes. We nevertheless think that 
our simulations have caught the essential properties of the topol- 
ogy of the epoch of reionization. Our simulations suggest that the 
idea that reionization proceeds strictly inside-out from beginning to 
nearly to the end may need revision if reionization indeed occurs in 
a photon- starved regime as suggested by the Lya forest data. 
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APPENDIX A: COMPARISON OF DIFFERENT 
METHODS OF GENERATING THE HALO FIELD 

In this Appendix, we compare three different ways of generating 
the density field and locating haloes within the simulation volume. 
Dark matter haloes were identified for density distributions with 
identical initial conditions within a simulation box of comoving 
length 50 Mpc with 256^ particles, giving a mass resolution of 
5.4 X 10^ Mq . For definiteness, we concentrate our comparison 
on z = 6 (which is the fiducial redshift of study throughout the 
paper). 

• N-body + FoF: In this approach, the dark matter density field 
is generated by running a full N-body simulation (with Gadgetll) 
and then a standard Friends-of-friends (FoF) algorithm with linking 
length h ^ 0.2 times the mean inter-particle separation is applied to 
find the haloes. Typically, one is able to identify haloes as small as 
^ 20 times the mass resolution which are consistent with theoreti- 
cal predictions of halo mass function. This is most accurate method 
to obtain the spatial distribution of dark matter haloes. The disad- 
vantage is that in order to achieve the dynamic range required for 
studying reionization is generally computationally expensive (both 
in terms of CPU time and memory). The density field and the loca- 
tion of dark matter haloes obtained in this way are shown in the top 
panels of Figure Al. 

• ZA + FoF: An alternate method of generating the density field 
is the Zel'dovich approximation. In this case, we have generated 
the density field at a given redshift by displacing the particles from 
their initial positions using the linear velocity field. This procedure 
is significantly less computationally expensive than a N-body sim- 
ulation and nevertheless gives give a reasonable representation of 
the density field at high redshifts. The location and mass of the 
haloes was then obtained with FoF halo finder with a variable link- 
ing length with h ^ 0.3 — 0.35. The detailed internal structure 
of the haloes is not correct in this case (the density profiles of the 
haloes is generally much more diffuse and the halo particle may 
even not not be bound). However, these details are not important 
for our work here where we want to investigate qualitatively the 
topology of reionization. The density field and the location of the 
haloes obtained in this way are shown in the middle panels of Fig- 
ure Al. One immediately appreciates that the visual impression of 
both the density structure and halo field generated by this approach 
is very similar to the previous one, the differences being rather mi- 
nor. 

• GRF + ES: The third method we have explored is evolving 
the initial Gaussian random field (GRF) linearly (i.e., multiplying 
by the appropriate growth factor) and applying the excursion set 
(ES) formalism to identify the haloes. The advantage in this case is 
that the formalism is computationally very cheap and can identify 
haloes as small as the mass resolution of the box. The disadvan- 
tage is that the linear density field does not necessarily capture the 
true density distribution which is a serious problem for the analyses 
presented here. The results obtained by this approach are shown in 
the bottom panels of Figure A 1. It is immediately apparent that the 
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density structure is drastically different from the previous two ap- 
proaches with no apparent filamentary networks visible. The same 
is true for the location of the haloes (though it should be mentioned 
that the number of haloes identified are much larger than the previ- 
ous methods as one can locate smaller haloes). A better match with 
the simulations can be achieved if both the densities and halo posi- 
tions are adjusted using the Zel'dovich approximation (Zahn et al. 
2007; Mesinger & Furlanetto 2007); however it is not clear how 
well the density peaks would correspond to halo locations if both 
are displaced independently. Since a reasonable representation of 
the density field and location of the haloes are vital for our work 
here, this very simple computationally least expensive scheme is 
unfortunately not appropriate for this work. 

The fact that the "ZA + FoF" method gives a reasonable approx- 
imation of the density and halo field can also be seen quantita- 
tively from Figure A2 where we have plotted the volume- weighted 



density distribution Pv(A) (left panel) and the power spectrum of 
density fluctuations P{k) (right panel) for the three methods. The 
density distribution obtained with the "ZA + FoF" method (dashed 
curve) closely resembles that obtained with "Nbody + FoF" (solid 
curve), which is quite different from the gaussian distribution (dot- 
dashed curve) obtained with the "GRF + ES" method. Similarly, 
the plots of the power spectrum shows that the "GRF + ES" method 
deviates from the "Nbody + FoF" at scales ~ 10/i~^ Mpc, while 
the "ZA -I- FoF" method is reasonable down to scales of a few 
Mpc. At smaller scales, the "Nbody + FoF" method generates more 
power than the other two cases due to a correct treatment of non- 
linearities. It appears thus fair to say that the "ZA -i- FoF" method 
is a good approximation for scales > lh~^ Mpc, which should be 
sufficient for generating the ionization maps in this work. 
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Figure A2. Left panel: the volume-averaged probability distribution of the density field Py (A) obtained from N-body simulations (solid curve), Zel'dovich 
approximation (dashed curve) and gaussian random field (dot-dashed curve) respectively. Right panel: the power spectrum of density fluctuations obtained by 
the same three methods. 
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Figure Bl. The effect of box size on the photon mean free path for models 
HR and IR-1.0. The solid (dot-dashed) and the dashed (dotted) curves rep- 
resent the results for the 100/i~^ Mpc and 200/i~^ Mpc box, respectively, 
for the IR-1.0 (HR) model. 



APPENDIX B: NUMERICAL CONVERGENCE 

In this appendix, we discuss the effects of limited box size and 
mass resolution on our results. For simplicity, we shall keep our 
discussion focussed on the models HR and IR-1.0. 

In order to study the effect of box size, we have run a simu- 
lation with a box of length 200 Mpc (comoving) with 2000^ 
particles, thus giving the same mass resolution as our fiducial box. 
We find that the effect on quantities like ionized fraction Qi{/S) and 
the distribution of 21cm brightness temperature Pm(A2i) is negli- 
gible for all models. The only significant effect of a larger box size 
concerns the evolution of the photon mean free path Amfp (which is 
shown in Figure Bl) and, to some extent, the 21cm power spectrum 

Aii(fc). 

For model IR-1.0, we find no significant effect of the limited 
box size on the shape or amplitude of A|i (/c) other that we are able 
to probe larger scales with a larger box size. The mean free path 
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Figure B2. The halo mass function at z = 6 for the high-resolution simu- 
lation. The points with errorbars show the results from our simulation; the 
vertical errors correspond to the statistical uncertainties while the horizontal 
errors denote the bin size. The solid curve is the theoretical mass function 
of Sheth & Tormen (2002), with the fitting function adopted from Jenkins 
et al. (2001). 



Amfp is not affected by the limited box size for scales smaller than 
the box as can be seen by comparing the solid and dashed curves in 
Figure Bl). However, with our fiducial box size of 100 Mpc, 
it is not possible to probe the IGM when the mass-averaged neutral 
fraction 1-xf < 0.01. If the box size is doubled to 200 Mpc, 
we are able to probe a much smaller neutral fraction 1 — xf^ < 
0.002. This confirms the result that larger boxes are essential when 
reionization enters its final stages. 

The requirement for larger box sizes is more apparent for 
model HR, where we find that the limited box size affects the value 
of Amfp for scales about half the box size (dotted and dot-dashed 
curves in Figure Bl). In fact, we find that a box size of as large as 
100 Mpc is only sufficient for neutral fractions l—xf- > 0.25. 
This is not surprising as the HR model tends to produce large ion- 
ized regions whose growth can be affected seriously with a limited 
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Figure B3. Effect of resolution on the photon mean free path for model IR- 
1.0. The solid curve represents the lower resolution simulation with 1000^ 
particles, while the dashed curve is for the higher resolution simulation with 
2000^ particles. The box size in both cases is lOO/i-i Mpc. 



fiducial simulation model IR-1.0. At the very late stages of reion- 
ization, however, the mean free path in the two cases is similar. In 
fact, at large ionized mass fraction, the only structures to remain 
neutral have intermediate densities, which should be equally well 
probed by the two simulations with different resolution. 



box size. We come to similar conclusions when studying the 21cm 
power spectrum. However the differences are not as statistically 
significant as the number of points which are neutral decreases dur- 
ing the late stages of reionization. 

Finally, we present the effect of numerical resolution on our 
analyses. For this purpose, we have run a simulation box of length 
100 Mpc (comoving) with 2000^ particles, which gives a mass 
resolution of Mpart = 9.02 x lQPh~^MQ. Applying the FoF 
method with adaptive linking length on this distribution, we are 
able to locate haloes as small as 9.02 x lO^/i^^M©, thus achiev- 
ing sensitivities corresponding to haloes able to cool via atomic 
transitions. The mass function of haloes at z = 6 for this high- 
resolution simulation is shown in Figure B2; we have also shown 
the corresponding theoretical mass function (Jenkins et al. 2001) 
for comparison. The halo mass function agrees now very well with 
the theoretical expectation for an even larger dynamic range. 

In Section 3 we have shown that the spatial distribution of 
sources of ionizing radiation have a huge effect on the ionization 
fields. Thus, it is naturally expected that the ionization maps would 
be very different for a high resolution box if we include all the low- 
mass sources. However, our main concern is to study the resolution 
effects for an identical source distribution is identical. Keeping that 
in mind, we include only sources with M > 10^ h~^MQ so that 
the source distribution is statistically identical to that in our fiducial 
box. For the high resolution box, we smooth the density field to a 
grid-size of 0.5/i~^ Mpc corresponding to 200^ grid points in the 
box. 

The main effect of the resolution enters into our results 
through the recombination rate. Since it is dependent on the lo- 
cal density, we find that the rate is higher when we include high 
resolution (i.e., high density) pixels in the analysis. We would thus 
expect, for example, that the mean free path is smaller in the high 
resolution simulation (even when the source distribution is statis- 
tically similar). That is indeed the case as is shown in Figure B3 
where we have compared the high resolution simulation with the 



